Bienvenid@s a la tercera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en el diseño de experimentos, test de hipótesis y regresión lineal. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Slides de las clases:
Enlaces a videos de las clases:
En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 3 son las siguientes:
# Manipulación de estructuras
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(tidyr)
# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)
##
## Attaching package: 'plotly'
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
##
## Attaching package: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combine
En este apartado se verán conceptos que será necesario revisar para responder a las siguientes preguntas.
En clases se han visto diferentes tipos de test de hipótesis para demostrar una proposición sobre algún parámetro. Uno de los test vistos en clases es el Z-Test, el cual su distribución del test estadístico bajo la hipótesis nula se puede aproximar a una Gaussina. Para la aplicación de este test, resaltan los siguientes puntos:
Para calcular la significancia estadística al igual que con otros métodos esta se debe calcular como:
Luego, dependiendo del objetivo del test tenemos las metodologías one-sample y two-sample. Utilizaremos One-Sample cuando nuestro objetivo es comparar la media de una muestra con la media de la población. El Z-score del One-Sample se define como:
\[Z-score_{One-Sample} = \dfrac{\bar x - \mu}{\dfrac{\sigma}{\sqrt n}}\] Donde \(\bar x\) es la media de la muestra, \(\mu\) es la media de la población, \(\sigma\) es la desviación estándar de la población y \(n\) es el tamaño de la muestra.
Por otro lado, se utiliza Two-Sample cuando queremos comparar la media de dos muestras. El Z-score de Two-Sample se define con la ecuación:
\[Z-score_{Two-Sample} = \dfrac{(\bar x_2 - \bar x_1) - (\mu_1 - \mu_2)}{\dfrac{\sigma_1}{\sqrt n_1}+\dfrac{\sigma_2}{\sqrt n_2}}\] Donde \((\bar x_2 - \bar x_1)\) es la diferencia de las medias de la muestra, \((\mu_1 - \mu_2)\) la diferencia de las medias de la población, \(\sigma_{1,2}\) la desviación estándar de la población y \(n_{1,2}\) el tamaño de las muestras.
En la práctica aparece la necesidad de testear múltiples hipótesis (por ejemplo en biología se pueden utilizar múltiples grupos de control o querer estudiar múltiples resultados de un mismo experimento), de esta forma la primera idea es testear individualmente cada una de las hipótesis, el problema de este enfoque es que la probabilidad de que se obtenga al menos un resultado significante crece rápidamente (con un nivel de significancia \(\alpha = 0.05\) y \(20\) test ya se alcanza una probabilidad de \(64\%\) de tener resultados significantes por azar).
Una forma de corregir los inconvenientes del método anterior es utilizar el método de Bonferroni correction quien propone cambiar \(\alpha\) por \(\alpha/m\) (donde \(m\) es la cantidad de test de hipotesis realizados), esto resulta que las probabilidades de rechazar por error se mantengan bajas. De esta forma los p-valores obtenidos en un test de hipótesis y al utilizar Bonferroni correction, quedan dados por el producto de un \(p-valor_{i}\) y la cantidad de test realizados: \(\text{p-valor}_{i}*m\).
El objetivo de esta pregunta es programar la potencia de un test de hipótesis y observar como se comportan las la hipótesis nula v/s la alternativa para un Z-test. Con el desarrollo de este ejercicio, podrán visualizar las diferentes partes que conforman a un test de hipótesis, identificar que es el p-valor y evidenciar como varia la potencia de un test one-sample y two-sample al variar \(\alpha\).
Para recordar; sabemos que en estadística el concepto de potencia viene dado por:
\[Power = 1 - \beta\]
Donde \(\beta\) es la probabilidad de obtener un error de tipo II. Con esto, la potencia estadística viene a representar la probabilidad de rechazar la hipótesis nula cuando esta es falsa. O sea, la potencia de una prueba es la probabilidad de encontrar un resultado positivo dado que este existe. Una de las formas de representar la potencia de un test es a través del siguiente gráfico:
Del gráfico, es posible visualizar que a medida que aumenta la diferencia en la media de la población, se obtienen mayores valores de potencia estadística.
Recordada que es la potencia de un test de hipótesis, a continuación, usted deberá programar una función que sea capaz de obtener la potencia de un Z-test one-sample y two-sample. Para esto por favor considere los siguientes puntos:
function(n1=NULL, sigma1=0.5,
n2=NULL,sigma2=0.5, mu.Ha=0 ,
mu.True=0, alfa=0.05)
De los argumentos, tendremos que: \(n1\) representa la cantidad de datos para la muestra 1, \(sigma1\) es la desviación estándar de la muestra 1, \(n2\) la cantidad de datos para la muestra 2, \(sigma2\) la desviación estándar para la muestra 2, \(mu.Ha\) el mu del test de hipótesis y \(mu.True\) la media de la población real. Notar que la presencia de una segunda muestra solo es para el caso two-sample, para el caso one-sample el argumento de entrada \(n2\) debería ser nulo.
Codificada la función realice los siguientes experimentos:
\[ n1=16, sigma1=16, mu.Ha=100 , mu.True=Variar, alfa=0.05 \] \[ n1=16, sigma1=16, mu.Ha=100 , mu.True= Variar, alfa=0.01 \] \[ n1=16, sigma1=16, mu.Ha=100 , mu.True= Variar, alfa=0.1 \]
Se le recomienda que la variación se realice a través de un
for y grafique las curvas dentro de un mismo gráfico para
observar potenciales diferencias entre ellas.
Diseñe un experimento one-sample y visualice cómo se comportan las distribuciones normales de la hipótesis nula y la hipótesis alternativa al variar \(\alpha\).
Diseñe un experimento Two-sample y visualice cómo se comportan las distribuciones normales de la hipótesis nula y la hipótesis alternativa al variar \(\alpha\).
Para el diseño de experimentos y/o comprobación de sus métodos puede serles útiles (no hay problema si decide utilizar los mismos ejemplos):
A su vez, en estas fuentes se incluyen las visualizaciones esperadas de las distribuciones normales para cada experimento. Estos experimentos principalmente deben consistir en ver cómo cambia la potencia al cambiar alguno de los parámetros (por ejemplo \(\alpha\)), y se espera que se entregue un comentario acorde.
Respuesta
power.z.test <- function(n1=NULL, sigma1=0.5,
n2=NULL, sigma2=0.5, mu.Ha=0 ,
mu.True=0, alfa=0.05){
# Si n2 es NULL, se realiza un Z-test de una muestra (one-sample)
if(is.null(n2)){
# Calcular el valor crítico de Z para el alfa dado (solo cola superior)
Z_crit = qnorm(1 - alfa)
# Cálculo del denominador para el Z-score
denom = sigma1 / sqrt(n1)
# X_bar es el valor crítico en el que se rechaza la hipótesis nula H0
X_bar = Z_crit * denom + mu.Ha
# Cálculo del Z-score
Z = (X_bar - mu.True) / denom
# Potencia del test
Power = 1 - pnorm(Z)
# Definir los límites mínimo y máximo para el gráfico de las distribuciones
min_lim = min(rnorm(1000, mean = mu.Ha, sd = denom)) -
round(min(rnorm(1000, mean = mu.Ha, sd = denom))) %% 10
max_lim = max(rnorm(1000, mean = mu.True, sd = denom)) +
round(max(rnorm(1000, mean = mu.True, sd = denom))) %% 10
# Crear el gráfico de las distribuciones
plot <- ggplot(data.frame(x = c(min_lim, max_lim)), aes(x)) +
stat_function(fun = dnorm, args = list(mean = mu.Ha, sd = denom),
col='red') +
stat_function(fun = dnorm, args = list(mean = mu.True, sd = denom),
col='blue') +
stat_function(fun = dnorm, args = list(mean = mu.True, sd = denom),
xlim = c(X_bar, max_lim), geom = "area", fill='red') +
geom_vline(xintercept = X_bar, linetype="dotted", size=1) +
annotate("label", x=X_bar, y=Inf, label="alpha", vjust=2) +
theme_minimal() +
ggtitle("H0 vs Ha") +
xlab(expression(bar(X))) + ylab("Densidad")
# Si n2 no es NULL, se realiza un Z-test de dos muestras (two-sample)
} else {
# Calcular el valor crítico de Z para el alfa dado (solo cola superior)
Z_crit = qnorm(1 - alfa)
# Cálculo del error estándar
denom = sqrt((sigma1^2 / n1) + (sigma2^2 / n2))
# X_bar es el valor crítico en el que se rechaza la hipótesis nula H0
X_bar = Z_crit * denom + mu.Ha
# Cálculo del Z-score
Z = (X_bar - mu.True) / denom
# Potencia del test
Power = 1 - pnorm(Z)
# Definir los límites mínimo y máximo para el gráfico de las distribuciones
min_lim = min(rnorm(1000, mean = mu.Ha, sd = denom)) -
round(min(rnorm(1000, mean = mu.Ha, sd = denom))) %% 10
max_lim = max(rnorm(1000, mean = mu.True, sd = denom)) +
round(max(rnorm(1000, mean = mu.True, sd = denom))) %% 10
# Crear el gráfico de las distribuciones
plot <- ggplot(data.frame(x = c(min_lim, max_lim)), aes(x)) +
stat_function(fun = dnorm, args = list(mean = mu.Ha, sd = denom),
col='red') +
stat_function(fun = dnorm, args = list(mean = mu.True, sd = denom),
col='blue') +
stat_function(fun = dnorm, args = list(mean = mu.True, sd = denom),
xlim = c(X_bar, max_lim), geom = "area", fill='red') +
geom_vline(xintercept = X_bar, linetype="dotted", size=1) +
annotate("label", x=X_bar, y=Inf, label="alpha", vjust=2) +
theme_minimal() +
ggtitle("H0 vs Ha") +
xlab(expression(bar(X))) + ylab("Densidad")
}
# Devolver el gráfico y la potencia como una lista
return(list(plot=plot, power=Power))
}
# Niveles de significancia (alfa) que se van a utilizar en los experimentos
alfas <- c(0.05, 0.01, 0.1)
# Experimento one-sample:
# Media poblacional verdadera
mu_true <- 95
# Media bajo la hipótesis alternativa
mu_ha <- 100
# Para cada valor de alfa, se calculan y grafican las distribuciones de H0 y Ha
for (alfa in alfas) {
result <- power.z.test(n1=16, sigma1=16, mu.Ha=mu_ha, mu.True=mu_true, alfa=alfa)
print(paste("Gráfico para alfa =", alfa, "en el test one-sample"))
print(result$plot)
}
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## [1] "Gráfico para alfa = 0.05 en el test one-sample"
## [1] "Gráfico para alfa = 0.01 en el test one-sample"
## [1] "Gráfico para alfa = 0.1 en el test one-sample"
# Experimento two-sample:
# Media de la hipótesis alternativa para las dos muestras
mu_ha_two_sample <- 100
# Diferencia entre las medias poblacionales verdaderas de las dos muestras
mu_true_two_sample <- 95
for (alfa in alfas) {
result <- power.z.test(n1=16, sigma1=16, n2=16, sigma2=16,
mu.Ha=mu_ha_two_sample, mu.True=mu_true_two_sample, alfa=alfa)
print(paste("Gráfico para alfa =", alfa, "en el test two-sample"))
print(result$plot)
}
## [1] "Gráfico para alfa = 0.05 en el test two-sample"
## [1] "Gráfico para alfa = 0.01 en el test two-sample"
## [1] "Gráfico para alfa = 0.1 en el test two-sample"
Esta pregunta tiene como objetivo comprender como funciona un test de hipótesis y como deberíamos abordar la realización de múltiples test de hipótesis con datos reales.
La pregunta deberá ser desarrollada utilizando el dataset
marketing_campaign.csv. Con esto, deberá programar un
Z-test, con el cual estudiará a través de experimentos el
Income de personas con los grados académicos
Graduation, Master y PhD. Para
realizar esto considere la elaboración de los siguientes puntos de forma
secuencial:
Graduation, Master y PhD por
separado.
Por ejemplo para el caso de Graduation pueden generar estructuras de la siguiente forma:
| ID | Graduation |
|---|---|
| 5524 | 58138 |
| 2174 | 46344 |
| 4141 | 71613 |
| 6182 | 26646 |
| 965 | 55635 |
| … | … |
Donde los valores en la fila de Graduation representan los sueldos de las diferentes personas que conforman el dataset. Un punto importante a considerar es que los datos para los diferentes grados académicos poseen diferentes numero de datos (no se asusten por esto).
Programar el método Z-test con la metodología one sample y two sample, obteniendo los p-valores a través de las alternativas one-sided y two-sided. Para el caso de one-sided, cree una función capaz de obtener la cola menor y mayor de la gaussiana.
El calculo de las diferentes alternativas para calcular los p-valores deberá ser un argumento de su función, donde señalando ‘menor’,‘mayor’ (para los casos one-sided) y ‘two-sided’ deberá obtener el valor pertinente para cada caso.
Genere una función que permita realizar solo múltiples test del
tipo two-sample y aplique bonferroni correction a los p-valores
obtenidos. Notar que los múltiples test deberá realizar la comparación
entre todos los elementos de entrada, por ejemplo si deseamos comparar
los ingresos de Graduation, Master y
PhD, se deberían comparar los ingresos de
Graduation v/s Master, Graduation
v/s PhD y Master y PhD
Codificada las funciones, realice los siguientes experimentos con su función de test de hipótesis:
Compruebe si la media de los ingresos para la variable
Graduation es similar a 52000. Señale formalmente este
experimento y obtenga los p-valores para las alternativas one-sided y
two-sided.
Compruebe si la diferencia entre los ingresos de las personas con
el grado académico Graduation es cercana a cero en relación
a la recibida por los Master y PhD. Para este
punto utilice la función que le permite realizar múltiples test del tipo
two-sample.
Para estos experimentos usted deberá escoger un \(\alpha\), y con éste comentar respecto a los p-valores obtenidos.
Para los diferentes experimentos considere que la desviación estandar
de la población para los diferentes income son los
siguientes:
\[\sigma_{Graduation} = 28180\] \[\sigma_{Master} = 20160\] \[\sigma_{PhD} = 20615\]
Respuesta:
df = read.csv('marketing_campaign.csv', sep='\t')
head(df)
income_master <- df$Income[df$Education == 'Master']
income_phd <- df$Income[df$Education == 'PhD']
income_graduation <- df$Income[df$Education == 'Graduation']
# Se quitan los valores nulos para evitar errores
income_graduation <- na.omit(income_graduation)
income_master <- na.omit(income_master)
income_phd <- na.omit(income_phd)
# Implementación de Z-test one-sided y two-sided
z_test <- function(data1=NULL, sigma1=0.5, data2=NULL, sigma2=0.5,
mu.Ha=0, test.type = c('menor', 'mayor', 'two-sided'),
verbose=TRUE){
# Verificación de argumentos
if(length(test.type) != 1 || !(test.type %in% c('menor', 'mayor', 'two-sided'))){
stop("Por favor, escoge un tipo de Test válido: 'menor', 'mayor' o 'two-sided'.")
}
# One-sample Z-test
if(is.null(data2)){
n1 <- length(data1)
mean_sample <- mean(data1)
Z_score <- (mean_sample - mu.Ha) / (sigma1 / sqrt(n1))
# Cálculo del p-valor
if(test.type == 'menor'){
p_value <- pnorm(Z_score)
} else if(test.type == 'mayor'){
p_value <- 1 - pnorm(Z_score)
} else if(test.type == 'two-sided'){
p_value <- 2 * (1 - pnorm(abs(Z_score)))
}
# Salida opcional
if(verbose){
cat("\tOne-sample Z-Test:\n\nData analizada:", deparse(substitute(data1)),
"\nZ =", Z_score, "P-value =", p_value, "\n\n")
}
return(p_value)
} else { # Two-sample Z-test
n1 <- length(data1)
n2 <- length(data2)
mean1 <- mean(data1)
mean2 <- mean(data2)
Z_score <- ((mean1 - mean2) - mu.Ha) / sqrt((sigma1^2 / n1) + (sigma2^2 / n2))
# Cálculo del p-valor
if(test.type == 'menor'){
p_value <- pnorm(Z_score)
} else if(test.type == 'mayor'){
p_value <- 1 - pnorm(Z_score)
} else if(test.type == 'two-sided'){
p_value <- 2 * (1 - pnorm(abs(Z_score)))
}
# Salida opcional
if(verbose){
cat("\tTwo-sample Z-Test:\n\nData analizada:", deparse(substitute(data1)),
"y", deparse(substitute(data2)), "\nZ =", Z_score,
"P-value =", p_value, "\n\n")
}
return(p_value)
}
}
z.test.multiple_testing <- function(list_data, list_sigma, test.type = 'two-sided', alfa = 0.05){
# Combinaciones de los grupos a comparar
combinations <- combn(names(list_data), 2)
p_values <- c()
results <- data.frame()
# Iterar sobre todas las combinaciones posibles de pares de grupos
for(i in 1:ncol(combinations)){
data1 <- list_data[[combinations[1, i]]]
data2 <- list_data[[combinations[2, i]]]
sigma1 <- list_sigma[[combinations[1, i]]]
sigma2 <- list_sigma[[combinations[2, i]]]
# Realiza el Z-test para cada par de grupos
p_value <- z_test(data1, sigma1, data2, sigma2, mu.Ha = 0, test.type = test.type, verbose = FALSE)
p_values <- c(p_values, p_value)
# Almacena los resultados en un data frame
results <- rbind(results, data.frame(Grupo1 = combinations[1, i], Grupo2 = combinations[2, i],
P_Value = p_value))
}
# Evaluar significancia
results$Significant <- results$P_Value < alfa
return(results)
}
# Primer experimento:
sigma_graduation <- 28180 # Desviación estándar poblacional
# Test two-sided
p_value_two_sided <- z_test(data1 = income_graduation, sigma1 = sigma_graduation,
mu.Ha = 52000, test.type = 'two-sided')
## One-sample Z-Test:
##
## Data analizada: income_graduation
## Z = 0.8539824 P-value = 0.3931147
# Test one-sided (mayor)
p_value_one_sided_mayor <- z_test(data1 = income_graduation, sigma1 = sigma_graduation,
mu.Ha = 52000, test.type = 'mayor')
## One-sample Z-Test:
##
## Data analizada: income_graduation
## Z = 0.8539824 P-value = 0.1965574
# Test one-sided (menor)
p_value_one_sided_menor <- z_test(data1 = income_graduation, sigma1 = sigma_graduation,
mu.Ha = 52000, test.type = 'menor')
## One-sample Z-Test:
##
## Data analizada: income_graduation
## Z = 0.8539824 P-value = 0.8034426
# Resultados
cat("P-valor (two-sided):", p_value_two_sided, "\n")
## P-valor (two-sided): 0.3931147
cat("P-valor (one-sided mayor):", p_value_one_sided_mayor, "\n")
## P-valor (one-sided mayor): 0.1965574
cat("P-valor (one-sided menor):", p_value_one_sided_menor, "\n")
## P-valor (one-sided menor): 0.8034426
# Segundo experimento:
list_data <- list(Graduation = income_graduation,
Master = income_master,
PhD = income_phd)
list_sigma <- list(Graduation = 28180,
Master = 20160,
PhD = 20615)
# Realiza los múltiples tests
results <- z.test.multiple_testing(list_data, list_sigma, test.type = 'two-sided', alfa = 0.05)
# Muestra los resultados
print(results)
## Grupo1 Grupo2 P_Value Significant
## 1 Graduation Master 0.883967005 FALSE
## 2 Graduation PhD 0.006691735 TRUE
## 3 Master PhD 0.022366590 TRUE
El objetivo de este problema es estudiar como realizar múltiples test
de hipótesis simultáneamente. Para esto en primer lugar se estudiara el
método “intuitivo”, donde veremos sus limitantes y se comparará con el
método llamado Bonferroni correction, posteriormente se
realizará un estudio practico con el dataset
ratones.csv.
Un investigador se ha colocado en contacto con ustedes señalándoles que realiza diariamente test de hipótesis entre las muestras que toma día a día en su laboratorio. Con esto, al investigador le urge saber si realizar multiples test de hipótesis sin una corrección podría afectar la toma de decisiones. Para comprobar esto, les solicita comprobar matemáticamente como se comporta la probabilidad de obtener al menos un resultado significativos al azar de sus experimentos diarios. Para esto, les señala que la la probabilidad de obtener un experimento por azar puede ser simulado a través de los casos exitosos de una binomial (valores mayores a cero), donde el numero de observaciones son la cantidad de experimentos (\(m\)) y la probabilidad queda dada por \(\alpha\) del test.
A continuación, se entregan unas indicaciones mas especificas para desarrollar la pregunta:
data <- read.csv("ratones.csv",sep= ";", stringsAsFactors = T)
head(data)
# p.value
data$p.value <- as.numeric(gsub(",", ".", as.character(data$p.value)))
# Count rows where p.value < 0.05
cat("Rejected values without Bonferroni:", sum(data$p.value < 0.05, na.rm = TRUE))
## Rejected values without Bonferroni: 8
# p.value.Bonferroni
data$p.value.Bonferroni <- as.numeric(gsub(",", ".", as.character(data$p.value.Bonferroni)))
# Count rows where p.value.Bonferroni < 0.05
cat("\nRejected values with Bonferroni:", sum(data$p.value.Bonferroni < 0.05, na.rm = TRUE))
##
## Rejected values with Bonferroni: 2
Notamos que sin la corrección de Bonferroni, 8 valores son rechazados mientras que la corrección solo rechaza 2. Esto nos asegura robustez frente a la posibilidad de falsos positivos. Sin embargo, se aumenta la probabilidad de falsos negativos ya que se rechaza menos la hipótesis nula. Por lo tanto, es coherente concluir que el método de Bonferroni es útil si se quieren evitar errores de tipo I, más puede ser arriesgado para evitar errores de tipo II.
Respuesta Aquí:
probEmpirica <- function(alpha,m){
n <- 10000 # Cantidad de veces que se va a repetir el experimento para estimar la probabilidad
res <-rbinom(n, m, alpha) #Resultados de los experimentos
prob <- length(res[res>0]) / n # Probabilidad empírica
return(prob)
}
# Real probability
probReal <- function(alpha, m) {
return(1 - (1 - alpha)^m)
}
# Different values for m
m_values <- 1:50
alpha <- 0.05
# Calculate both probabilities
empirical_probs <- sapply(m_values, function(m) probEmpirica(alpha, m))
real_probs <- sapply(m_values, function(m) probReal(alpha, m))
# Graficar resultados
plot(m_values, empirical_probs, type = "o", col = "blue", pch = 16, ylim = c(0, 1),
xlab = "Número de ensayos (m)", ylab = "Probabilidad",
main = "Comparación de Probabilidad Empírica y Real")
lines(m_values, real_probs, type = "o", col = "red", pch = 17)
legend("bottomright", legend = c("Probabilidad Empírica", "Probabilidad Real"),
col = c("blue", "red"), pch = c(16, 17), lty = 1)
Notamos que los resultados son muy parecidos, realizando 10000
experimentos para la probabilidad empírica. Si se disminuye el número de
experimento es cierto que aumenta la variabilidad de esta, pero sigue la
misma tendencia de ir convergiendo a 1 a medida que aumenta el número de
ensayos.
El valor de alpha solo influye en la velocidad con la que ambas probabilidades convergen a 1, pero independientemente de este valor, siempre se produce que ambas probabilidades convergen eventualmente a 1 a la misma velocidad. Si se disminuye alpha convergen más lento y si se aumenta convergen más rápido. Esto nos indica que si se podrían realizar más experimentos para obtener resultados significativos, pero el precio a pagar sería reducir el nivel de significancia.
No es útil, ya que uno no quiere tener una alta probabilidad de obtener resultados significativos por azar, y este método claramente expone que al aumentar los ensayos aumenta esta probabilidad (Y mucho).
# Código modificado utilizando Bonferroni correction
probEmpiricaBonferroni <- function(alpha,m){
actual_alpha <- alpha / m
n <- 10000 # Cantidad de veces que se va a repetir el experimento para estimar la probabilidad
res <-rbinom(n, m, actual_alpha) #Resultados de los experimentos
prob <- length(res[res>0]) / n # Probabilidad empírica
return(prob)
}
# Real probability
probRealBonferroni <- function(alpha, m) {
actual_alpha <- alpha / m
return(1 - (1 - actual_alpha)^m)
}
# Different values for m
m_values_b <- 1:50
alpha_b <- 0.05
# Calculate both probabilities
empirical_probs_b <- sapply(m_values_b, function(m) probEmpiricaBonferroni(alpha_b, m))
real_probs_b <- sapply(m_values_b, function(m) probRealBonferroni(alpha_b, m))
# Graficar resultados
plot(m_values_b, empirical_probs_b, type = "o", col = "blue", pch = 16, ylim = c(0, 1),
xlab = "Número de ensayos (m)", ylab = "Probabilidad",
main = "Comparación de Probabilidad Empírica y Real con Bonferroni Correction")
lines(m_values_b, real_probs_b, type = "o", col = "red", pch = 17)
legend("bottomright", legend = c("Probabilidad Empírica", "Probabilidad Real"),
col = c("blue", "red"), pch = c(16, 17), lty = 1)
Vemos que la probabilidad se mantiene consistenmente cercana a 0 a
medida que aumenta m, lo cual nos indica que el método de Bonferroni si
es efectivo para evitar la obtención de falsos positivos. El problema es
que si se aumenta mucho m, aumenta la probabilidad de obtener errores de
tipo II, es decir, falsos negativos, ya que el método puede ser muy
conservador (Solicitando un valor de significancia innecesariamente
bajo).
El objetivo de la siguiente pregunta es aplicar los conceptos de regresión lineal vistos en clases para implementar desde cero un función capaz de realizar una regresión simple y múltiple.
Para este problema, ustedes deberán estudiar el comportamiento de los
clientes de un holding de salud. Para esto, se les hace entrega del
dataset insurance.csv para que estudien la creación de un
modelo lineal con sus datos. Antes de comenzar a trabajar, se señalan
las diferentes variables que componen al dataset:
Es importante que considere que cada una de las filas representa un cliente distinto para el holding.
Dentro del estudio, el holding de salud le solicita estudiar los
comportamientos de los clientes fumadores y no fumadores, por lo que se
le aconseja separar el dataframe original en fumadores y no fumadores.
En el estudio, realicen un modelo lineal que tiene como variable de
respuesta a charges y los datos que mejor se correlacionan
para los clientes fumadores y no fumadores. Para esto, deberán realizar
las siguientes actividades:
charges.lm(), ya que esto le servirá para observar
si la regresión lineal es significativa.Nota: No esta permitido utilizar comandos que obtengan los valores solicitados directamente a menos que se le permita en la pregunta.
Las variables escogidas son bmi para fumadores (ya que es un indicador de salud) y age para no fumadores (ya que los cargos tienden a aumentar con la edad, independiente de si fuman o no).
# Parte I.
# Separate smokers from non-smokers
insurance <- read.csv("insurance.csv", stringsAsFactors = T)
smokers <- insurance[insurance$smoker == "yes",]
non_smokers <- insurance[insurance$smoker == "no",]
# Calculate cors
cor_smokers <- cor(smokers %>% select(age, bmi, children, charges))
cor_smokers
## age bmi children charges
## age 1.00000000 0.05967388 0.08118329 0.36822444
## bmi 0.05967388 1.00000000 -0.01261916 0.80648061
## children 0.08118329 -0.01261916 1.00000000 0.03594501
## charges 0.36822444 0.80648061 0.03594501 1.00000000
cor_no_smokers <- cor(non_smokers %>% select(age, bmi, children, charges))
cor_no_smokers
## age bmi children charges
## age 1.00000000 0.12263798 0.03339533 0.62794678
## bmi 0.12263798 1.00000000 0.01920841 0.08403654
## children 0.03339533 0.01920841 1.00000000 0.13892870
## charges 0.62794678 0.08403654 0.13892870 1.00000000
# Linear model
lm_reg <- function(X, Y) {
n <- length(Y)
# Calcular las sumas necesarias
sum_X <- sum(X)
sum_Y <- sum(Y)
sum_XY <- sum(X * Y)
sum_X2 <- sum(X^2)
# Calcular los coeficientes beta_1 (pendiente) y beta_0 (intersección)
beta_1 <- (n * sum_XY - sum_X * sum_Y) / (n * sum_X2 - sum_X^2)
beta_0 <- (sum_Y - beta_1 * sum_X) / n
# Devolver los coeficientes
return(c(beta_0, beta_1))
}
# FFunction to make predictions
predict_lm <- function(X, coefficients) {
return(coefficients[1] + coefficients[2] * X)
}
# Adjust model
coefficients_smokers <- lm_reg(smokers$bmi, smokers$charges)
cat("\nIntercept (beta_0) smokers:", coefficients_smokers[1], "\n")
##
## Intercept (beta_0) smokers: -13186.58
cat("Slope (beta_1) smokers:", coefficients_smokers[2], "\n")
## Slope (beta_1) smokers: 1473.106
coefficients_non_smokers <- lm_reg(non_smokers$age, non_smokers$charges)
cat("\nIntercept (beta_0) smokers:", coefficients_non_smokers[1], "\n")
##
## Intercept (beta_0) smokers: -2091.421
cat("Slope (beta_1) smokers:", coefficients_non_smokers[2], "\n")
## Slope (beta_1) smokers: 267.2489
get_R_squared <- function(SSE, SST, Y, X, coefs) {
preds <- predict_lm(X, coefs)
SSE <- sum((Y - preds)^2)
SST <- sum((Y - mean(Y))^2)
r_squared <- 1 - (SSE / SST)
n <- length(Y)
p <- 1
r_squared_adjusted <- 1 - ((1 - r_squared) * (n - 1) / (n - p - 1))
return(c(r_squared, r_squared_adjusted))
}
smokers_r_squared <- get_R_squared(SSE_smokers, SST_smokers, smokers$charges, smokers$bmi, coefficients_smokers)
non_smokers_r_squared <- get_R_squared(SSE_non_smokers, SST_non_smokers, non_smokers$charges, non_smokers$age, coefficients_non_smokers)
cat("\nR² smokers:", smokers_r_squared[1], "\n")
##
## R² smokers: 0.650411
cat("R² adjusted smokers:", smokers_r_squared[2], "\n")
## R² adjusted smokers: 0.6491257
cat("\nR² non smokers:", non_smokers_r_squared[1], "\n")
##
## R² non smokers: 0.3943172
cat("R² adjusted non smokers:", non_smokers_r_squared[2], "\n")
## R² adjusted non smokers: 0.3937468
# Smokers scatterplot
ggplot(smokers, aes(x = bmi, y = charges)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, color = "blue") +
labs(title = "Regresión lineal: Gastos Médicos vs. IMC (fumadores)",
x = "Índice de Masa Corporal (BMI)",
y = "Gastos Médicos (Charges)") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
# Non smokers scatterplot
ggplot(non_smokers, aes(x = age, y = charges)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE, color = "blue") +
labs(title = "Regresión lineal: Gastos Médicos vs. Edad (No fumadores)",
x = "Edad (Age)",
y = "Gastos Médicos (Charges)") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
# Parte II
# Usamos la regresión multivariable con representación matricial
# Selección de las variables independientes (predictoras) y dependiente
X_smokers <- smokers %>% select(age, bmi)
Y_smokers <- smokers$charges
X_non_smokers <- non_smokers %>% select(age, bmi)
Y_non_smokers <- non_smokers$charges
# Función para calcular la regresión lineal múltiple
lm_reg_multi <- function(X, Y) {
X <- as.matrix(cbind(Intercept = 1, X))
Y <- as.matrix(Y)
# Calcular los coeficientes usando la fórmula de mínimos cuadrados
beta <- solve(t(X) %*% X) %*% t(X) %*% Y
return(beta)
}
# Entrenar el modelo para fumadores
coefficients_smokers <- lm_reg_multi(X_smokers, Y_smokers)
cat("\nCoeficientes para fumadores:\n")
##
## Coeficientes para fumadores:
cat("Intercepto (beta_0):", coefficients_smokers[1], "\n")
## Intercepto (beta_0): -22367.45
cat("Edad (beta_1):", coefficients_smokers[2], "\n")
## Edad (beta_1): 266.2922
cat("BMI (beta_2):", coefficients_smokers[3], "\n")
## BMI (beta_2): 1438.091
# Entrenar el modelo para no fumadores
coefficients_non_smokers <- lm_reg_multi(X_non_smokers, Y_non_smokers)
cat("\nCoeficientes para no fumadores:\n")
##
## Coeficientes para no fumadores:
cat("Intercepto (beta_0):", coefficients_non_smokers[1], "\n")
## Intercepto (beta_0): -2293.632
cat("Edad (beta_1):", coefficients_non_smokers[2], "\n")
## Edad (beta_1): 266.8766
cat("BMI (beta_2):", coefficients_non_smokers[3], "\n")
## BMI (beta_2): 7.075477
# Función para hacer predicciones con múltiples variables
predict_lm_multi <- function(X, coefficients) {
X <- as.matrix(cbind(Intercept = 1, X))
return(X %*% coefficients)
}
# Función para calcular R² y R² ajustado
get_R_squared_multi <- function(Y, X, coefs) {
preds <- predict_lm_multi(X, coefs)
SSE <- sum((Y - preds)^2)
SST <- sum((Y - mean(Y))^2)
r_squared <- 1 - (SSE / SST)
n <- length(Y)
p <- ncol(X)
r_squared_adjusted <- 1 - ((1 - r_squared) * (n - 1) / (n - p - 1))
return(c(r_squared, r_squared_adjusted))
}
# Evaluar el modelo para fumadores
r_squared_smokers <- get_R_squared_multi(Y_smokers, X_smokers, coefficients_smokers)
cat("\nR² para fumadores:", r_squared_smokers[1], "\n")
##
## R² para fumadores: 0.7532403
cat("R² ajustado para fumadores:", r_squared_smokers[2], "\n")
## R² ajustado para fumadores: 0.7514192
# Evaluar el modelo para no fumadores
r_squared_non_smokers <- get_R_squared_multi(Y_non_smokers, X_non_smokers, coefficients_non_smokers)
cat("\nR² para no fumadores:", r_squared_non_smokers[1], "\n")
##
## R² para no fumadores: 0.3943673
cat("R² ajustado para no fumadores:", r_squared_non_smokers[2], "\n")
## R² ajustado para no fumadores: 0.3932257
summary(lm(smokers$charges ~ smokers$bmi + smokers$age))
##
## Call:
## lm(formula = smokers$charges ~ smokers$bmi + smokers$age)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14604.4 -4315.1 -240.5 3638.0 29316.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -22367.45 1931.86 -11.58 <2e-16 ***
## smokers$bmi 1438.09 55.22 26.05 <2e-16 ***
## smokers$age 266.29 25.06 10.63 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5754 on 271 degrees of freedom
## Multiple R-squared: 0.7532, Adjusted R-squared: 0.7514
## F-statistic: 413.6 on 2 and 271 DF, p-value: < 2.2e-16
summary(lm(non_smokers$charges ~ non_smokers$bmi + non_smokers$age))
##
## Call:
## lm(formula = non_smokers$charges ~ non_smokers$bmi + non_smokers$age)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3237.7 -1947.5 -1354.7 -650.6 24430.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2293.632 804.061 -2.853 0.00442 **
## non_smokers$bmi 7.075 23.877 0.296 0.76703
## non_smokers$age 266.877 10.245 26.048 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4669 on 1061 degrees of freedom
## Multiple R-squared: 0.3944, Adjusted R-squared: 0.3932
## F-statistic: 345.4 on 2 and 1061 DF, p-value: < 2.2e-16
Notamos que aumenta considerablemente el valor de R^2 llegando a 0.75 en el caso de la regresión multivaribale. Esto se debe a la agregación de una nueva variable que aporta mayor información del comportamiento. Algo a notar si es que para el caso de no fumadores el valor de R^2 no varía, lo que nos indica que no existe mayor correlación entre la variable cargos y las variables escogidas juntas.
Sí sería recomendable utilizar los modelos para la predicción ya que estos coinciden con los entregado por el comando lm() y la evaluación es mejor en términos de la métrica de R^2.
A work by CC6104